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Abstract 

We present an explicit solution based on the phase-amplitude approximation of the Fokker-Planck equa- 
tion associated with the Langevin equation of the birhythmic modified van der Pol system. The solution 
enables us to derive probability distributions analytically as well as the activation energies associated to 
switching between the coexisting different attractors that characterize the birhythmic system. Comparing 
analytical and numerical results we find good agreement when the frequencies of both attractors are equal, 
while the predictions of the analytic estimates deteriorate when the two frequencies depart. Under the ef- 
fect of noise the two states that characterize the birhythmic system can merge, inasmuch as the parameter 
plane of the birhythmic solutions is found to shrink when the noise intensity increases. The solution of the 
Fokker-Planck equation shows that in the birhythmic region, the two attractors are characterized by very 
different probabilities of finding the system in such a state. The probability becomes comparable only for a 
narrow range of the control parameters, thus the two limit cycles have properties in close analogy with the 
thermodynamic phases. 
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The van der Pol oscillator is a model of self-oscillating system that exhibits periodic os- 
cillations. A modified version - essentially a higher order polynomial dissipation - has 
been proposed as a model equation for enzyme dynamics. This model is very interesting 
as a paradigm for birhythmicity, it contains multiple stable attractors with different natural 
frequencies, therefore it can describe spontaneous switching from one attractor to another 
under the influence of noise. The noise induced transitions between different attractors de- 
pend upon the different stability properties of the attractors, and are usually investigated by 
means of extensive Langevin simulations. We show that the associated Fokker-Planck equa- 
tion, in the phase-amplitude approximation, is analytically solvable. The phase amplitude 
approximation requires a single frequency, and therefore fails when the two frequencies of 
the birhythmic system are significantly different. However, the approximation is not severe, 
for it explains the main features of the system when compared to the numerical simulations 
of the full model. The approximated Fokker-Planck equation reveals the underlining struc- 
ture of an effective potential that separates the different attractors with different frequency, 
thus explaining the remarkable differences of the stability between the coexisting attractors 
that give rise to birhythmicity. Moreover, it reveals that the noise can induce the stochastic 
suppression of the bifurcation that leads to birhythmicity. Finally, the approximated solu- 
tion shows that the system is located with overwhelming probability in one attractor, thus 
being the dominant attractor. Which attractor is dominant depends upon the external con- 
trol parameters. This is in agreement with the general expectation that in bistable systems 
the passage from an attractor to the other resembles phase transitions, since only in a very 
narrow interval of the external parameters it occurs in both directions with comparable 
probabilities. 



I. INTRODUCTION 



A stochastic dynamical system is a dynamical system under the effects of noise. Such effects 
of fluctuations have been of interest for over a century since the celebrated work of Einstein yj]. 
Fluctuations are classically referred to as "noisy" or "stochastic" when their suspected origin im- 
plicates the action of a very large number of variables or degrees of freedom. For a linear system 
this leads to the phenomenon of diffusion, while the coupling of noise to nonlinear deterministic 
equations can lead to non-trivial effects |2Ll3j]. For example, noise can stabilize unstable equilibria 
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and shift bifurcations, i.e. the parameter value at which the dynamics changes qualitatively |4j,|5|]. 
Noise can also lead to transitions between coexisting deterministic stable states or attractors such 
as in birhythmic or bistable system [60 . Moreover, noise can induce new stable states that have no 
deterministic counterpart, for instance noise excites internal modes of oscillation, and it can even 
enhance the response of a nonlinear system to external signals QVL LSD - 

In this paper, we investigate analytically the effects of an additive noise on a special bistable 
system that displays birhythmicity - coexisting attractors that are characterized by different fre- 
quencies I2I-I6D. We examine a birhythmic self-sustained system described by the modified van- 
der Pol oscillator, subjected to an additive Gaussian white noise 116 1 . 

Our main aim is to use the phase-amplitude approximation II 1711 . a standard technique for van 



der Pol linn and van der Pol - like systems II 1811 . to derive an effective Fokker-Planck equation 



11911 that can be analytically managed. This allows us to analytically derive the activation energies 
associated to the switching between different attractors !(], I^]]. The analytical solution of the 
approximated model is not limited to vanishingly small noise intensity as it was done for the 
numerical estimate of the escape time [6J] to derive the pseudopotential [21]. Another purpose of 
the present paper is to verify, with numerical simulations, that in spite of the approximations the 
analytical probability distribution is reliable. 

The paper is organized as follows: Section II presents the modified van-der Pol system with 
an additive Gaussian white noise. Section III deals with the derivation and analysis of an effec- 
tive Fokker-Planck equation for the birhythmic modified van der Pol oscillator. The probability 
distribution given by the approximated Fokker-Planck equation is analyzed and the activation en- 
ergies are derived. In Section IV, we integrate numerically the stochastic second order differential 
equation and discuss the results. Section V concludes. 



H. THE BIRHYTHMIC PROPERTIES OF THE NOISY MODEL 
A. The modified van der Pol oscillator with an additive noise 

The model considered is a van der Pol oscillator with a nonlinear dissipation of higher poly- 
nomial order described by the equation (overdots as usual stand for the derivative with respect to 
time) 

x — — x 2 + ax 4 — (3x 6 )x + x = 0. (1) 
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This model was proposed by Kaiser Ii22l1 as more appropriate than the van der Pol oscillator to 
describe certain specific processes in biophysical systems. In fact the modified van der Pol-like 
oscillator described by Eq. (OQ) is used to model coherent oscillations in biological systems, such 
as an enzymatic substrate reaction with ferroelectric behavior in brain waves models (see Refs. 
|23|426|] for more details). From the standpoint of nonlinear dynamics, it represents a model which 
exhibits an extremely rich bifurcation behavior. The quantities a and (3 are positive parameters 
which measure the degree of tendency of the system to a ferroelectric instability compared to 
its electric resistance, while ji is the parameter that tunes nonlinearity [23]. The model Eq.© 
is a nonlinear self-sustained oscillator which possesses more than one stable limit-cycle solution 
Ii27ll . Such systems are of interest especially in biology, for example to describe the coexistence 
of two stable oscillatory states, as in enzyme reactions 112 811 . Another example is the explanation 
of the existence of multiple frequency and intensity windows in the reaction of biological systems 



23, 



29- 



3211 . Moreover, the 



when they are irradiated with very weak electromagnetic fields H24I . 
model under consideration offers general aspects concerning the behavior of nonlinear dynamical 
systems. Kaiser and Eichwald |32j] have analyzed the super-harmonic resonance structure, while 
Eichwald and Kaiser 12211 have found symmetry-breaking crisis and intermittency. 

In Ref. Il23n an analytical approximation has been derived for the coexisting oscillations of the 
two attractors with different natural frequencies for the deterministic part of the model equation. 
A numerical investigation of the escape times (and hence of the activation energies) has suggested 
that the stability properties of the attractors can be very different Q6Q . It has further been shown 
that time delayed feedback leads to stabilization [330, also in the presence of external noise B20Q. 

Noise can enter the system for instance, through the electric field applied to the excited enzymes 
which depends on the external chemical influences or through the flow of enzyme molecules. 
One can therefore assume that the environmental influence contains a random perturbation and to 
postulate that the activated enzymes are subject to a random excitation governed by the Langevin 
version of Eq. (OQ), namely: 

x- -x 2 + ax 4 - /3x 6 )x + x = T(t). (2) 

T(t) can be assumed to be an additive Gaussian white noise with arbitrary amplitude D and 
it has the properties: 



< r(t) >= 

< T(t),T(t') >= 2D5(t-t') 



(3) 



which completely determine its statistical features. The noise term is here treated as external 
Ii36ll . i.e. due to a disturbance from the environment and not subject to the fluctuation dissipation 
theorem. 



B. Birhythmic properties 

Without noise (T = 0), Eq.© reduces to the modified version of the van der Pol oscillator 
© which has steady- state solutions that depend on the parameters a, (3 and [i and correspond to 
attractors in state space. The dynamical attractors of the free-noise modified van der Pol Eq.© 
have been determined analytically, the expressions of the amplitudes Aj and frequency flj (i= 1,2,3) 



2311 . in which the periodic solutions 



of the limit-cycle solutions have been established in Ref. 
of the modified van der Pol oscillator © are approximated by 

x(t) = A cos tit. (4) 

The amplitude A is independent of the coefficient \i up to corrections of the order /j 2 and implicitly 
given by the relation: 

^A e - -A* + -A 2 - 1 = 0. (5) 

64 8 4 

The coefficient fi enters in the expression for the frequency Cl as a second order correction: 

Q — 1 + fi 2 uj2 + o(/i 3 ) (6) 
thus the deviations of the frequency from the linear harmonic solution are characterized by an 



amplitude dependent frequency Q26Q - 



93/3 2 , 19 69aB . in ,67/3 3a 2 . ,« ,73/3 a . . fi ,1 a . 3 <9 

u 2 = — —A -A 10 + ( — — H )A 8 -( — — H )A 6 + ( 1 )A A A 2 (7) 

65536 16384 v 8192 1024 y v 2048 96 ; v 128 24 y 64 

Depending on the values of the parameters a and (3, the modified van der Pol oscillator possesses 
one or three limit cycles. In fact, Eq.Q can give rise to one or three positive real roots that 
correspond to one stable limit cycle or three limit cycle solutions (of which two are stable and one 
is unstable), respectively. The dynamical attractors and birhythmicity (i.e. the coexistence between 
two stable regimes of limit cycle oscillations) are numerically found solving the amplitude Eq.© 
If]]. The three roots Ai, A 2 , and A 3 denote the inner stable orbit, the unstable orbit, and the outer 
stable orbit, respectively. When three limit cycles are obtained, Eq.© supplies the frequencies 



^1,2,3 in correspondence of the roots Ai^. Being one of the attractors unstable, the system only 
displays two frequencies Qi ;3 (and hence birhythmicity) at two different amplitudes Ai t3 , while 
the unstable limit cycle of amplitude A 2 represents the separatrix between the basins of attraction 



of the stable limit cycles. We s 



low in Fig.l the region of existence of birhythmicity in the two 



parameters phase space (a- (3) 11231. 12611 (the two coexisting stable limit cycle attractors can be 



found in [6J]). The question we want to address is the influence of noise on the above properties 
investigating the response of an additive Gaussian white noise in the phase-amplitude limit. In 
Ref.fl] the system has been numerically tackled in the regime of vanishingly small noise. In this 
limit the escape rate gives an effective potential that acts as an activation barrier. We employ the 
phase-amplitude approximation that should be both faster (being analytical) and more accurate at 
finite values of the noise, as will be discussed in the next Sect. III. 



III. ANALYTICAL ESTIMATES 

The analytic results on the deterministic system are based on the approximated cycle given by 
Eq.©. Quite naturally, one can treat the noise in the system starting from such approximation. 
To this extent, we rewrite the Langevin Eq.© in a system of two coupled first order differential 
equations: 

x = u, 

u = n(l - x 2 + ax 4 - (3x 6 )u -x + T. (8) 
We seek for solution in the context of the phase-amplitude approximation, i.e. letting the amplitude 



and the phase of Eq.© to be time dependent 



lli: 



x = A(t) cos(ttt + <f>(t)) 

u = -A(t)uj sm(Qt + (j)(t)). (9) 

Inserting Eq.© into Eq.®, one retrieves a system of two Langevin equations for the amplitude 
A(t) and phase variables, that is, of course, as difficult to manage as the original model ©. 



We will follow the standard analysis of nonlinear oscillators j 18L 13411 that consists in assuming that 



in a period 2n/{l the variables A(t) and (pit) do not change significantly, so one can average the 



effect of the random perturbation Il20ll . Although in principle this method also relies on the small- 



ness of the noise, since the averaging requires that the approximate solution © is not significantly 
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altered in a cycle 2-k/VL, the procedure has proven very robust in a similar van der Pol - Duffing 



oscillator It 1 811 . It is important to note that for fj, — the system reduces to the harmonic oscillator, 
as described by the solution Eq.© with constant amplitude and phase. Since we are interested in 
the influence of noise D and nonlinear dissipation (a and (3) in birhythmic systems, we keep the 
parameter /i small (/i = 0.1). If the present model is employed to model the population of enzyme 
molecules, the parameter represents the difference between the thermal activated polarization and 



the external field induced polarization 112711 . However, we note that for a birhythmic system a fur- 
ther difficulty occurs: the system has two different frequencies ^ f2 3 , while the approximation 
© is monorhythmical. Assuming that the two frequencies are not too different, we insert Eq.© 
into Eq.® and average, to retrieve the effective (and simpler) Langevin equation for the amplitude 
A and phase <p variables: 

. q 2 - 1 r~D~ _ . 

' = "XT - V 2A& < m >' 



A = ^[(i-^ + ^-^)]-^<r(t)>. do) 



We thus study the system in the slow averaged variables; for the slow variables the average noise 



can still be considered white and uncorrelated I117h . and the Fokker-Planck equation associated to 
the Langevin model (PTOl) reads: 

dP dS^ dS A 



dt ' d(f) OA' 
where S = + S A is the probability current defined by: 

= KfP - -^-(Kp^P), 



(11) 



S A = KfP-^K^P). (12) 
The drift coefficients K± and K A associated to Eq.dTQl) read: 

Kt - 

K A = ^[1 - - A 2 + —aA 4 - -f3A 6 } + (13) 
1 2 1 4 8 64^ J 2Q 2 A v J 

The off diagonal diffusion coefficients Kp A and K A ^ vanish, while the diagonal coefficient read: 



L2 (fL4) 2 ' 



We seek for stationary solutions, dP/dt = of Eq.dTTb. We note that in the averaged equation 
© for the phase A the phase (j) does not appear, and therefore the integration over all phases gives 
rises to a normalization constant. We therefore only seek solutions for the probability distribution 
associated to the constant probability current, S A = const. Moreover, since the probability distri- 
bution must vanish for A = oo, we can set the constant to 0. Finally, the equation for the radial 
part of the probability distribution P reads: 



or, explicitly: 



S A^ Q ^ K A p= d_ {K A,A pl (15) 



P(A) = cAe W {^[l - l -A* + ±aA* - A/L4 6 ]}, (16) 



where c is a constant of normalization. This solution contains as particular cases the harmonic 
oscillator (p, < 0, a = j3 — 0, and discarding the A 2 / 8 term) and the standard van der Pol 
oscillator (/i > 0, a — f3 — 0). 

The probability distribution is in general very asymmetric, for most of the parameters a or (3 
one can localize the probability function around a single orbit. Before proceeding further in our 
analysis, it should be noted that the peaks of the probability distribution can be located using the 
following equation: 



dlog(P) =Q 



A2 ~^-°- (17) 



dA 

For D — 0, the amplitude (flTT) coincides with the deterministic amplitude equation (5) ||20fl . In 
Fig. 2 we report the influence of the noise intensity D on the region of multi-limit cycle orbits of 
Fig. 1. In the parametric (a, /3)-plane of Fig. 2 it is evident the effect of the noise intensity D on the 
transition boundary between the appearance of single and multi-limit cycles orbits: the bifurcation 
that leads to birhythmicity is postponed under the influence of noise [35]. As a consequence, the 
region of existence of three limit cycles, a condition for birhythmicity, decreases with the increase 
of the noise intensity and disappears altogether for high noise intensity. 

An important feature of birhythmicity in the present model is highlighted in Fig. 3. We first 
define 7\ 3 of the probability to find the system in the basin of attraction of each stable orbit 1 and 
3: 

-Ai 



Vi = [ * P(A)dA, 
Jo 



POO 

V 3 = / P{A)dA. (18) 
Ja 2 



These quantities measure the relative stabilities pertaining to the attractors 1 and 3 and are related 
to the resident time by the relation Vi, 3 = T^/fTi + T 3 ), so that T>i/P 3 = T^jT^. In Figfjwe 
show in the parameter plane a — (3 the locus where the system stays with equal probability on both 
attractors (solid line) T\ = T 3 . We also show two further curves: the limit where the first attractor 
is much more stable than the other (T x > 10T 3 , circles) and the passage to the reverse situation 
(T 3 > lOXi, crosses). From the figure it is evident that the outer attractor is dominantly visited in 
most of the parameter plane. Moreover, the transition from the two opposite cases (i.e., a change 
of two order of magnitudes of the relative resident times) occurs with a very narrow change of the 
control parameters a and (3. The drastic change is further investigated in Fig. 4, where we show 
a blow-up of the crossover region around 7\ = T 3 for different values of the parameter (3. The a 
value is increased up to the maximum value when birhythmicity disappears. The general behavior 
observed for all (3 values, closely reminds phase transitions: the probability to find the system in 
one condition (around the attractor A{) or the other (around the attractor A 3 ) drastically changes 
in a very small interval of the a parameter. The same behavior, this time with a constant value of 
a and varying (3 is shown in Fig. 5. The effect of the noise intensity is much less pronounced, 
see Fig. 6. It is apparent that the effective temperature is capable to cause a crossover between 
the residence times only in a very narrow region of the phase space, inasmuch as the noise causes 
a contraction of the region of existence of birhythmicity. However it is evident that the transition 
is much slower, and a crossover only occurs in the limited parameter space around a = 0.05, 
f3 = 0.0005. 

The stability properties of the two attractors have also been investigated in the limit of small 
noise values [6], where it has been found the same asymmetrical behavior of the probability dis- 
tribution, with a sudden change for small variations of a and (3. In fact one can notice that the 
effective Langevin equation (flOl) amounts to the Brownian motion of a particle in a double well, 



whose potential reads [20J]: 



F A (A) = -^[(1 - -A 2 + —aA 4 - —(3A% (19) 
v 1 4 [V 8 24 256 n 

It is therefore evident that the transition from the inner orbit A\ to the outer orbit A 3 through 
the unstable orbit A 2 , as well as the inverse process, can be interpreted as the diffusion over an 
effective potential barrier, and therefore the escape times are given by the Kramer's inverse rate 
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1 36], for instance used in Josephson physics to detect the classical quantum/classical crossover 
0711 or for signal detection 113 811 . 

4ft 2 



ti^ 3 oc exp 
r 3 ^i oc exp 



D 

4ft 2 



(F A (A 2 ) - F A {A X )) 
(F A (A 2 ) - F A (A 3 )) 





'AUi' 


= exp 


D 




~AU S 


= exp 


D 



(20) 



Thus the average time to pass from one attractor to the other is analogous to the passage over 
a barrier. The pseudopotential barrier numerically derived in Ref.|^] is therefore, in the phase- 
amplitude approximation, an effective potential for the amplitude variables [2(3]. Since the effec- 
tive potential is analytical, we can confirm several features of the pseudo-potential, for instance 
that the potential barriers are proportional to the nonlinear parameter [i [6!]. It is also interest- 
ing to investigate the behavior of the potential barriers of Eq.(l20l) as a function of the parameters 
a an 0, the analogous of the analysis of Eq.(fT6l) in Figs. 3,4,5,6. Inspection of the effective 
potential (1201 ), confirms that it is very asymmetrical, since one energy barrier is generally much 
higher than the other. Combining this observation with the exponential behavior of the escape 
rates (l20l) one deduces that the system does not equally stays on both attractors, but rather it 
clearly exhibits a preference for one attractor with respect to the other (the relative occupancies 
read T 1 /T 3 ~ exp[(A{7 3 — AU\)/D] Q36D). One concludes that the birhythmic system behaves 
as a bistable tunnel diode yfjfl: keeping fixed a control parameter (say /3) and changing the other 
(a in this case) the weight of the probability distribution is concentrated in the proximity of one 
or the other of the two stable deterministic solutions of Eq. (OQ), thus obtaining again a first order 
phase transition. This result supports the notion that the analogy with phase transitions is generic 
for bistable oscillators n39fl . 



IV. NUMERICAL SIMULATIONS AND RESULTS 



To check the validity of the approximations behind the analytic treatment that has led to the so- 
lution (fT6l) . we have performed numerical simulations of the Langevin dynamics ©■ There are 
several methods and algorithms for solving second-order stochastic differential equations ikOll as 
the implicit midpoint rule with Heun and Leapfrog methods or faster numerical algorithms such as 
the stochastic version of the Runge-Kutta methods and a quasisymplectic algorithm ||4l|l. To prove 



that the simple procedure given by the Euler algorithm is reliable, we have employed it in a 



:ew 



selected points with two different methods. The starting point is the Box-Mueller algorithm I142ll 
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to generate a Gaussian white noise distributed random variable T^t from two random numbers a 
and b which are uniformly distributed on the unit interval [0, 1]. The random number approximates 
the effect of the noise of intensity D over the interval At in the Euler algorithm for the integration 
of Eq.®. We have then halved the step size At until the results became independent of the step 
size; the step size used for all numerical integration is At = 0.001. To verify the numerical re- 



sults obtained with the Euler method, we have used a quasi- symplectic algorithm of Mannella 114 lh 



to numerically compute the probability distribution. The logic behind the choice to compare the 
Euler algorithm with a quasi-symplectic algorithm is that the nonlinear dissipation of the model 
CO) oscillates and vanishes twice in each cycle. We have therefore checked the results with an 

o 

algorithm that has proved to perform independently of the dissipation value 114011 . 

In Fig. 7 we plot the behavior of the probability distribution P as a function of the amplitude 
A for several values of the noise intensity D, when the frequencies of both attractors are similar, 
i.e. fi x ~ fi 3 ~ 1. It clearly shows that the system is more likely found at two distinct distances 
from the origin, the essential feature of birhythmicity. In general, for the set of parameters a = 
0.083, B = 0.0014, the probability distribution P is asymmetric. As observed in Sec. Ill the 

nn 

probability distribution changes with a small variation of the parameters a and 6 [a, I20J]. It is 
important to note that the agreement between numerical and analytical results is fairly good for 
low A values around the inner orbit, when the frequency of one attractor is very similar to 1, 
while for larger amplitude A the agreement becomes progressively poorer. However, it seems that 
the phase-amplitude approximation is capable to capture the main feature of the phenomenon: 
an increase or a decrease of the amplitude when the fluctuation parameter D is varied. At high 
fluctuations (D > 1) the system becomes monorhythmical, see also Fig. 2, thus confirming the 
noise induced transition from bimodal to unimodal, sometimes referred to as phenomenological 
bifurcations II 1 8M . 

As mentioned in Sec. Ill, the phase-amplitude approximation is not appropriate when the two 
frequencies of the attractors are different, i.e. Vt\ ^ f2 3 . In fact numerical simulations in this case 
show a poor agreement, see Fig. 8 where VL\ ~ 1 and Q 3 ~ 0.8. This shows the limitations of this 
analysis of phase- amplitude approximation. 

Let us return to Eq.(fT7l) that shows how the orbits radii depend on the noise intensity D. The 
analytical and numerical behaviors of the limit cycle attractors are reported in Figs. 9 and 10 
that show amplitudes A\ 3 and the associated bandwidths AA 1)3 (the width when the height of the 
probability peaks is reduced of a factor 2) as a function of the noise intensity D for two sets of 
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parameters a and 0. We find that the amplitudes A\ and A 3 change very slightly when the noise 
intensity increases. Also the bandwidth slightly increases with the noise intensity D. Through Eqs. 
( 1161171 ) one can derive the behavior of the effective potential barriers iB I^oll . We have numerically 
compared the analytic predictions with simulations in Fig. 11, where we plot AUi and A?7 3 as a 
function of the parameter a with (3 = 0.002. It should be noted that according to the numerical 
results of Ref . D6Q , varying a from a = 0.095 to a = 0.135 the system passes from the region 
where Qi ~ f2 3 to the region with Q\ ^ f2 3 . It is clear that in general the two energy barriers 
are very different. For low a values A?7 3 is well approximated by the analytic approach, while 
for larger a the agreement becomes progressively poorer. Nevertheless it seems that the phase- 
amplitude approximation is capable to capture the main feature of the phenomenon: an increase 
or a decrease of the activation energies when the dissipation parameters are varied. An analogous 
behavior is observed in Fig. 12, where we plot the behaviors of AU\ and AU3 as a function of the 
parameter (5. 



V. CONCLUSIONS 



We have approached a theoretical description of the temporal evolution of the modified van der 
Pol oscillator with an additive Gaussian white noise in the region where birhythmicity (in the ab- 
sence of noise) occurs. To get an analytical insight on this system we have used an explicit solution 
based on the phase-amplitude approximation of the Fokker-Planck equation to analytically derive 
the probability distributions. The activation energies associated to the switches between different 
attractors have been derived analytically and numerically. We have found that the agreement is 
fairly good. The characteristics of the birhythmic properties in a modified van der Pol oscillator 
are strongly influenced by both the nonlinear coefficients a, (3 and the noise intensity D. The 
boundary of the existence of multi-limit-cycle solutions, in the parametric (a, /3)-plane, decreases 
with the increase of the noise intensity D. Finally, the analytic estimate of the stability of the two 
attractors varies with the control parameters (the dissipation a and /3) in a way that resembles phase 
transitions: for most parameters value the system is located around only one attractor, the other 
being visited with a vanishingly small probability. Only at special values of the control parameters 
the residence times are comparable, in agreement with experimental observations of birhythmicity 
in Biological systems: the passage from an attractor to another only occurs by varying the external 



parameters and not under the influence of noise [|43L 
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FIG. 1: Parameters domain for the existence of a single limit cycle (white area) and 
three limit cycles (black area) solutions ofEq. (Q/or /! = 0.1. 
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FIG. 2: Effect of the noise intensity D on the boundary between the region of one 
and three limit-cycle solutions in the parametric (a, j3)-plane of the Fokker-Planck 
Eq.[T7\) for \i = 0.1 as in Fig\j\ 
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FIG. 3 : Behavior of the residence times in the parameter space. The solid line denotes 
the locus Ti = T 3 , while circles and crosses denote the situation where T\ = 10T 3 
and T\ = T3/IO, respectively. The dashed line is the border of existence ofbirhyth- 
micity. The noise level is D = 0.1 
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FIG. 4: Residence tunes as a function of the parameter a for different values of the parameter (3. The noise 
level is D = 0.1,, the nonlinearity fx = 0.1 
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FIG. 5: Residence times as a function of the parameter ft for different values of the parameter a. The noise 
level is D = 0.1, the nonlinearity fi = 0.1 
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FIG. 6: Residence tunes as a function of the noise intensity D for different values of the parameter a. The 
second dissipation parameter reads /3 = 0.0005, the nonlinearity fx = 0.1. 
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FIG. 7: Asymmetric probability distributions for different values of the noise intensity D versus the ampli- 
tude A when the frequencies of both attractors are identical i.e fii ~ ~ 1. Parameters of the system are 
fj, = 0.1 and a = 0.083, /3 = 0.0014. 
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FIG. 8: Probability distribution versus the amplitude A when the frequencies of the attractors are not 
identical i.e 0i ^ 0,$. Parameters of the system are D = 0.1, yU = 0.1, (i): a = 0.09, (3 = 0.0012, Qi ~ 
l,fl 3 ~ 0.85 and(ii): a = 0.1, /3 = 0.014, n x ~ 1,0 3 - 0.8. 
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FIG. 9: Variation of the amplitudes At and the bandwidths AAi versus the noise intensity D. Lines and 
symbols denote analytical and numerical results, respectively. The circles and dot-dashed lines refer to the 
inner attractor A\, solid lines and triangles to the outer attractor A3. The parameters used are fx = 0.1, 
a = 0.1, f3 = 0.002. 
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FIG. 10: Variation of the amplitudes Ai and the bandwidths AAi versus the noise intensity D. Lines and 
symbols denote analytical and numerical results, respectively. The circles and dot-dashed lines refer to the 
inner attractor A\, solid lines and triangles to the outer attractor A3. The parameters used are fi = 0.1, 
a = 0.12, /3 = 0.003. 
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FIG. 12: Behavior of energy barriers versus /3. Solid lines denote the analytical results, while dashed lines 
with triangles denote numerical simulations. Parameters of the system are \i = 0.1 and a = 0.13. 
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